Exact and approximate methods of calculating the sum of states for noninteracting 
classical and quantum particles occupying a finite number of modes 
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We present exact expressions for the sum of states of noninteracting classical and quantum parti- 
cles occupying a finite number of modes with arbitrary spacings. Exploiting a probabilistic analogy, 
we derive an analytic fourth-order approximation to the density of states, which captures its vari- 
ance and kurtosis, and is superior to the previous, commonly used methods for all three particle 
statistics. Our approach employs a simple exact method of calculating the moments of the micro- 
canonical density of states for quantum particles, which requires less computational effort than the 
commonly used saddle-point approximation. We test our methods numerically and discuss their 
applicability to various physical systems. 



I. INTRODUCTION 



Erwin Schrodinger wrote "There is, essentially, only 
one problem in statistical thermodynamics: the distribu- 
tion of a given amount of energy E over TV identical sys- 
tems" [1]. As the number of particles in a system grows 
to infinity, almost all its copies (or, expressing the same 
notion differently, almost all its configurations) have the 
same internal energy. In this limit, the canonical ensem- 
ble (where we control the temperature) is equivalent to 
the microcanonical ensemble (where we control the en- 
ergy). The opportunity to pass between these different 
descriptions gives the possibility to choose the smarter 
way to measure a given quantity. Very often the calcula- 
tion of thermodynamical quantities in the microcanonical 
ensemble is an impractical task and thus one is forced to 
resort to the canonical ensemble. However, their equiva- 
lence breaks down for small or perfectly isolated systems. 
One then has to work with the microcanonical ensemble, 
whose workhorse is the sum of states function, which 
counts the number of microstates realizing a particular 
value of the control parameter. This is the case, for ex- 
ample, when investigating atoms in microcavities [3 or 
elongated magneto-optical traps Q , heavy nuclei [J Q , 
coupled spins on a lattice Q or calculating the entropy 
of a black hole 0,1- 

In all physical problems listed above the system is 
composed of N classical (distinguishable), bosonic or 
fermionic particles which can occupy S + 1 modes, num- 
bered from to S. Each mode can be (^-degenerate and 
has the excitation E s > 0, with E s < E s+ i and Eq = 0. 
The last assumption does not lead to the loss of gener- 
ality, because we can always shift the excitations up or 
down to ensure that the lowest mode has zero excitation. 
In many applications, we have E s = s (e.g. lattice spins 
or harmonic traps) or E s = s 2 (square potential wells). 

The control parameter for which we calculate the sum 



of states is defined as 



S gs-l 
s=0 t=0 



(1) 
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where n s t = 0, 1, . . . is the number of particles occupy- 
ing mode s and (in the case of degenerate modes) its 
submode t, subject to the constraint 

S Ss-l 

J2J2n st =N. 

s=0 t—0 

The number of physical states of N particles realizing a 
given M is called the sum of states tl(N, M, S). When 
particles are distinguishable, M can be the total magne- 
tization of a system of lattice spins @ or the total area 
of a black hole [7]. For bosons, M is often the total en- 
ergy of particles in a harmonic oscillator potential 0, H| 
or atoms in a microwave cavity 0. For fermions, it can 
be the z component of the total angular momentum of a 
heavy nucleus [J, [5[ . For the sake of generality, we call 
M a total excitation. 

The paper concerns with the calculation of tl(N, M, S) 
for noninteracting classical, bosonic or fermionic particles 
occupying a finite number of modes with arbitrary spac- 
ings. First, we investigate exact (the popular recursive 
and proposed iterative) expressions for the marginal or 
cumulative sum of states, which are useful in practical 
calculations. Exploiting a probabilistic analogy, we then 
obtain an analytic fourth-order approximation to the nor- 
malized sum of states (density of states) u>, which cap- 
tures its mean, variance and kurtosis. The approximation 
is extended to an analytic-numerical scheme which can 
also match the skewness of to and is better suited to han- 
dle excitation spectra with strongly non-homogeneous 
densities. We demonstrate numerically that our analytic 
approximation is superior to the previous ones (the Gaus- 
sian one and its polynomial corrections fitted also to the 
kurtosis — e.g. P. M liol fllj ) and matches closely the exact 
sum of states for all three particle statistics. Additionally, 
our approach employs a simple exact method of calculat- 
ing the moments of the microcanonical density of states 



2 



for quantum particles, which requires less computational 
effort than the commonly used saddle-point approxima- 
tion (see e.g. d, [Hj]). We discuss the applications of our 
methods to different physical systems pointing out their 
advantages over commonly employed approaches. 



II. EXACT EXPRESSIONS FOR fl 

We begin by reviewing the known recursive exact ex- 
pressions for the sum of states of classical and quantum 
noninteracting particles 1 1 .'ii . and proceed later to derive 
their iterative counterparts. 



A. Recursive expressions 

1 . Sum of states 

The total excitation of N particles distributed among 
5 + 1 modes is defined by (|T|). If the mode 5 (degenerate 
or not) is occupied by n particles, than the remaining 
excitation M—nEg can be redistributed among the N—n 
particles in w(N, n)Cl(N—n, M—nEs, 5—1) ways, where 
the factor 



w(N, n) 



N\ 
(N-n)ln\ 



(classical) 
(quantum) 



accounts for the particle statistics. On the other hand, if 
the mode 5 is gs-degenerate, the n particles occupying 
it can be reshuffled in 



{gs) n (classical) 
r+flT 1 ) ( b °sonic) 
X 9 n) l ^<gs (fermionic) 



ways [14j |. Hence, we have the recursive expression 

Sl(N, M, 5) = 

■n(gs,N) 

w(N,n)v(g s ,n)to(N - n,M - nE s ,S - 1) 

n=0 

where 



(2) 



N (classical/bosonic) 
min(gs,N) (fermionic) 



The boundary condition for Q is 

n(N,M,0) = 5 Mfi v(g ,N). 



(3) 



For many values of n the factor Vt(N — n, M—nEs, S — 
1) will be zero, because n will be either too high or too 
low. Thus, it makes sense to derive tighter bounds for 
the range of summation over n. For its given value M, 



the number n of particles occupying mode s cannot be 
higher than 

Qs = min([M/E s \ ,v(gs,N)) 

(where [x\ is the highest integer number less than or 
equal to x), otherwise the total excitation would exceed 
M. 

On the other hand, for classical particles or bosons 
the maximum M that N — n particles occupying modes 
below S can realize is (N — ri)Eg_±, which leads to the 
inequality 

M < nE s + {N-n)E S -i = NE S -i+n(E s -E S -i) ■ (4) 
Since Es > Es-i, n is greater or equal to Ps defined as 
~M-NE S - 



Ps = max 



Es-i 



, 0] (class. /bosonic) (5) 



(6) 



where |~x] is the lowest integer number greater than or 
equal to x. Equation ([2]) is modified to 

ft(N, M, S) = 

Qs 

J2 w{N, n)v(g s , n)Q(N - n, M - nE s , S - 1) . 

n=P s 

The remaining question is the definition of Ps for 
fermions, for which it should be higher due to the Pauli 
exclusion principle. In principle, we can derive an in- 
equality similar to ((U), but it would be tedious to use, 
as it would depend on the excitation values of multiple 
modes s < S. However, we can safely use the Ps as 
defined for bosons and classical particles ([5]) when deal- 
ing with fermions, at the cost of summing over a larger 
number of zero values. 

For classical particles, we can also use recursion in the 
number of particles, not in the number of modes, writing 



Q(N, M, S)=£ g s Q{N -1,M-E„,S).. 



(7) 



s=0 



with the same boundary condition as before. This for- 
mula holds an advantage over ((SJ), because it does not 
require the evaluation of factorials. 



2. Cumulative sum of states 

We define the cumulative sum of states as the number 
of states with total excitation less than or equal to M, 
S(iV, M, S). To distinguish between the two quantities, 
the previously defined sum of states Vt(N, M, S) can also 
be called the marginal sum of states. S(iV, M, S) satis- 
fies a recursive equation similar to ©, but without the 
lower limit Ps (which ensured that the total excitation 
constraint ([!]) was satisfied), 



E(N, M, S) 



71 = 



w(N, n)v(g s ,n)E(N -n,M - nE s , 5-1) 
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with a boundary condition also somewhat different 
from 

X(N,M,0) = v(g ,N). (8) 

Analogously to ([7]), we obtain for classical particles 
only 

s 

£(AT, M, S) = J2 9sZ{N - 1, M - E S ,S) , 

s=0 

with the same boundary condition 

In numerical computations, the recursive formulas ([6]) 
and J7| require the storage of the intermediate values of 
17 or S. This means that the excitation values E s and M 
must be discretized, which for irregular excitation values 
leads to large memory requirements. On the other hand, 
they do not iterate over every configuration of the system, 
avoiding the problem of the exponential growth of their 
total number. 



B. Iterative expressions 



2. Cumulative sum of states 

The difference between Q(N, M, S) and E(N, M, S) is 
the relaxation of constraint Q when calculating the lat- 
ter. Hence, we can easily write an iterative exact formula 
for it, similar to ©, 

Qs Qi 
E(N,M,S)= E E u ({^})- 

ns — ni— 

The other conditions and definitions remain unchanged. 

III. ANALYTIC APPROXIMATIONS FOR ft 

Formulas for the sum of states , © an d ^ , though 
accurate, are not very convenient for analytic calcula- 
tions. We will now derive two analytic approximations 
for 17 (iV, M, S) (more precisely, for its version ui(N, M, S) 
normalized to unity), by exploiting an analogy between 
a system of classical/quantum independent particles and 
a set of independent /correlated random variables. 



Sum of states 



A. Classical particles 



In certain applications, it is necessary to have access 
to each of the system configurations counted by the sum 
of states. For this purpose, the recursive formulas can be 
converted to an explicit iterative expression for the sum 
of states: 

Qs Qi 

n(N,M,s)= E ••• E U (W)> ( 9 ) 

ns = Ps «l=fl 

where n s = J^I'Zo n ss', 



P, = max 



E s — E s -i 



,0 



M-El s+1 E t n t 



E, 



and n = N - J2t=i n t- 
The factor 



,v(9s,N 



E *)) . 

!=s+l / 



S / S \ 

o({n s }) = JJ w I N - ^ n s ,n s \v(g s 

s=0 \ t=s+l ) 



n s ) 



is the number of particle states realizing a given combi- 
nation {n s }. Equation (jH]) provides a simple way to sum 
over all combinations of {n s }, subject to the constraints 

J2s=o Us = N an< ^ E s =o n sE s = M, and takes into ac- 
count the mode degeneracies. Again, the difference be- 
tween the formulas for bosons and fermions amounts to 
excluding, in the latter case, all combinations where at 
least one n s > g s . 



When dealing with classical particles, one can associate 
the excitation state of each particle with an independent, 
uniformly distributed random variable Mj taking real 
values Eq,..., Es, each having a probability 1/ X^=o 9 s > 
so that its ex 
and variance 



so that its expected value E[M,-] = Yls=o 9sE s / 2f=o 9 



Var[Mj] = E[(Mj - E[M,]) 2 ] = 



■■- _ Es=o9s(E s -E[M 3 ]r 



The total excitation M corresponds to the sum of these 
variables, M = YljLi Mj and the probability that M — 
M is given by P(M = M) = u(N, M, S). The variance 
of M is thus 



JV 



Var[M] =^Var[Mj-] 



and the mean 



N 



E[M] = 5>[M,-] . 

3=1 



(10) 



(11) 



In the case of equal mode spacings, E s = s, and g s = 1 
we have 

E[M] = ^ and Var[M] = NS ^ + l \ 

The numerical examples presented below correspond to 
this case, but the presented method will be applicable to 
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any reasonably homogeneous excitation spectrum (i.e. in 
which the mode spacings do not systematically increase 
or decrease over the range of modes considered, leading 
to the skewness of the sum of states function). (We will 
discuss the non- homogeneous spectra in Sec. IIII CI ) 

From the Central Limit Theorem it follows that the 
probability distribution of M is well approximated for 
large TV by a Gaussian distribution $^0-2 with mean 
given by (jlll) and variance a 2 by (fit))) . Hence, 



n(N, M, 5) 
(5+1) 



1 



N 



27rVar[M] 



: exp 



(M — E[Af]) 2 
2Var[Af] 



which can be also expressed as 



fl(N, M, 5) « n(N, E[M],S) exp - 



(M - E[M] 



2 Vax[M] / 

(12) 

In Fig. [T] we compare the exact, obtained using 
the methods from the previous section, and approxi- 
mate results for Q(N,M,S), plotting the logarithm of 
ui(N, M, 5) to highlight better the discrepancies for ex- 
treme values of M. There is an excellent agreement 
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FIG. 1: Logarithms of exact and approximate values of 
oj(N, M, S) for N = 50 and S = 5 (classical particles, equal 
mode spacings). 

between the Gaussian approximation and exact results 
around M = NS/2. However, further away from the 
mean their logarithms differ, suggesting that the higher 
moments of the exact density of states are different from 
the Gaussian one. Indeed, the excess kurtosis of a single 
particle variable Mj is equal (for E s = s and g s = 1) to 



KurtfMj 



hence 



Var[Sj] 2 



Kurt[M] = —Kurt^] 



12 + 65(5' + 2) 
55(5 + 2) 

12 + 65(5 + 2) 
5iV5(5 + 2) 



(13) 



The exact probability distribution of M as a random vari- 
able is platykurtic (has negative excess kurtosis, i.e. thin 
tails), as opposed to the Gaussian distribution, which 
has zero kurtosis. For large N, Kurt [M] —> 0, which 
is consistent with its distribution approaching Gaussian. 
For smaller N, such as 50, the negative excess kurtosis 
is still non-negligible and needs to be accounted for. A 
much better fit to the exact sum of states is obtained by 
adding a (M— E[M]) 4 term to the exponent of (|12p. yield- 
ing Cl(N,M, 5) ~ exp(-a(Af-E[Af]) 2 -6(M-E[M]) 4 ). 
A numerical least-squares fit of this form to exact data 
is displayed in Fig. [T] as "4th order (fitted)" curve. 

To derive the analytic approximation including fourth- 
order terms, we study the properties of the following 
probability density 



*p(z) 



2V2a 



;V32a ff iir l/4 (l/32a) 



exp 



z 

2^2 



ax 



(14) 



where K n (x) is the modified Bessel function of the 
second kind. Our model for oj(N,M,S) is therefore 
u(N, M, S) = <p(M — E[M]). The case a = corresponds 
to the Gaussian density. For N — 50 and 5 = 5, the 
best-fit value of a (displayed in Fig. [TJ is 0.004933, thus 
we will concentrate our attention on the a —¥ limit, and 
derive the expressions for a and a which replicate best 
the shape of the tu(N, M, 5) function for each N and 5, 
a(N 7 5) and a(N, 5). The mean of the distribution 
is zero, its variance 



v(a,a) w cr 2 (l - 12a) 



(15) 



and the the fourth moment about the mean 

.4 , _ 4, 



/Lt 4 (cr, a) 



V(z)z^dz w cr 4 (3 - 96a) 



Hence, the excess kurtosis is a function of a only 

, , m{a,a) „ 24a(l + 18a) 
K[a) — — 777 — 3 — 



v(a, a) 



(1 - 12a) 



For a given value of Kurt[M] < 0, we can solve for a 
such that (fLT| has the same excess kurtosis, obtaining 



a(N, S) 



Kurt[M] 



1 - 5 Kurt [M] - 1 



12(3 + Kurt[M]) 
Inserting this into P^)) . we obtain 



Var[M](3 + Kurt[Af]) 
cr(iV,5)= =■ , 

\ 4- ^/l-5Kurt[Af] 

and the approximation reads 

Q(M, N, 5) Q(N, E[M] , 5) x 

/ (M — E[Af]) 2 , AT m (A/-E[M])^ 



(16) 



(17) 



(18) 



2a(N, Sy 



a(7V,5) 4 



5 



The variance and kurtosis of M are given by (fTU|) 
and (|T5|) . As shown by the "4th order (analytic)" curve 
in Fig. [TJ for N — 50 the analytic approximations to 
a(N, S) and a(N, S) work better than the Gaussian ap- 
proximation, but are not optimal. However, one can ex- 
pect that their accuracy will grow with N due to de- 
creasing kurtosis. To test this, we plot the exact values 
and both (Gaussian and fourth-order) analytic approxi- 
mations for N = 1000 and S = 5 in Fig. [5] Indeed, the 
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FIG. 2: Logarithms of exact and approximate values of 
aj(iV, M, S) for N — 1000 and S = 5 (classical particles, equal 
mode spacings). The analytic fourth-order formula is much 
closer to the exact calculations than in the case of TV = 50 
(Fig- ED- 

fourth-order analytic approximation works much better 
than the Gaussian and matches quite closely the exact 
calculations. 



B. Quantum particles 

Given the values of the mean, variance and kurtosis 
of M, we can use the Gaussian or fourth-order approxi- 
mation for N independent fermionic or bosonic particles. 
However, the computation of these quantities is not as 
easy as for classical particles, because the particles' ex- 
citations cannot be identified with independent random 
variables, due to their quantum nature. For example, 
given a system of two bosons which can occupy states 
and 1, the joint probability distribution of their exci- 
tations in the microcanonical ensemble is given by the 
table 



Configuration Probability 



(0,0) 
(1,0) or (0,1) 
(1,1) 



1/3 
1/3 
1/3 



Such a probability distribution cannot be realized by two 
independent random variables. On the other hand, the 



mode occupation numbers n s t are independent, if we as- 
sume that the number of particle in the system is not 
fixed. Let us assume that we have two modes or 1 and 
that their occupation numbers n\^n% are independent 
and uniformly distributed between and 2. There are 9 
combinations of occupation numbers, each having prob- 
ability 1/9. If we now impose the constraint n\ +ri2 = 2, 
we obtain a conditional distribution of occupation num- 
bers given in the table below (where we have used the 
curly braces to distinguish particle configurations from 
mode occupations): 



Occupations Probability 



{2,0} 

{1,1} 
{0,2} 



1/3 
1/3 
1/3 



The distribution in the second table is identical to the 
one in the first table, with {2,0} corresponding to con- 
figuration (0,0), {1,1} to configurations (1,0) or (0,1), 
and {0,2} to (1,1). It follows that the sum of states 
Q.{M, N, S) for quantum particles is proportional to the 
probability distribution of the variable M\N — N (M 
conditional on N equal to N), where M is the total exci- 
tation of the system (J^ at E s n s t), and N (also a random 
variable now) is the total number of particles (J2st n st)- 
The fact that the independent random variables are now 
associated with mode occupation numbers and not with 
particle excitations is a reflection of the fact that the 
particles themselves are excitations of a matter field in 
Quantum Field Theory, while the constraint on the num- 
ber of this excitations (N — N) is a conservation law 
characterizing massive particles. 

To approximate the distribution of M\N — N ana- 
lytically, we calculate its mean, variance and kurtosis, 
and then use the fourth-order approximation (jT5J . The 
calculation of the above moments is done recursively. 
Let M t denote the total excitation of first t quantum 
modes. (An excitation mode E s with degeneracy g s is 
equivalent to g s quantum modes, each with the same 
excitation value.) Let Hk(n,t) = K[Mj^\N = n] and 
p(n, t) = P(N = n) conditioned on particles occupying 
first t quantum modes. We have 



N 



p(n, t) ~ Yl P( n ~ n '> 1 ~ > H P( n > = 1 ( 19 ) 



n'=0 



n=0 



and 



/i fe (n,f)= ^E[(£ t m + M t _i) fc |iV = n-m 

m— 

x p(n — m, t — 1) . 



(20) 



For example, 

n 

/ii(n, t) = ^2 [Etm + /ii(n — to, t — l)]p(n — m,t — 1) 



m=0 
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FIG. 3: Logarithms of exact and approximate sum of states 
for bosons in a ID harmonic potential, with N — 100 and 150 
modes. 



In general, we can calculate p(n,t) from p(n,t — 1) and 
Hk(n, t) from fj,^(n,t — 1). The starting conditions are 
/ife(n,0) = and p(n,0) = 5% t o- Unlike in the case of 
classical particles, the above approximation requires a 
recursive numerical calculation. However, this calcula- 
tion is less onerous than the one in ([2]), as we do not 
loop over M values. Additionally, our exact method of 
calculating the moments fXk(n,t) is much simpler than 
the commonly used non-exact approach, which requires 
the approximation of the microcanonical density of states 
using the canonical or grand canonical density of states 
and saddle-point approximation [l2j ■ 

We test the approximation in the case of a ID har- 
monic oscillator potential, for which g s = 1 and E s — s 
(equal mode spacings). In the case of bosons, the approx- 
imation agrees well with exact results (Fig. [3]). Similar- 
results were obtained for fermions. Figure|3]compares the 
results of the approximation for different particle statis- 
tics, showing how important it is to include the effect 
of quantum indistinguishability and especially the Pauli 
principle (which is omitted by e.g. the Bethe formula used 
to describe the sum of states of nuclear spins [3, [Bj]). In 
the tests for non-equal mode spacings the proposed ap- 
proximation is at least as good as the Gaussian one. 



C. Non- homogeneous excitation spectra 

The formula (|18|) describes a symmetric sum of states 
function. Thus, our approximation will not work too well 
for systems with strongly non-homogeneous excitation 
spectra, where the dependence of the sum of states on 
M is strongly skewed, e.g. with E s — s n for n > 2. This 
can be improved at the cost of making the approximate 




M 



FIG. 4: Approximate sum of states for bosons, fermions and 
classical particles in a ID harmonic potential, with N = 100 
and 150 modes. 



density of states non-analytic, and using the distribution 

[exp (- - si^ail) Z<ZQ 
VskcwM- ) iz J of a2(z l A (21) 

[ eX P(, 2^1 ^ ) Z>Z 

instead of (|14[) . Together, the parameters zo, a± t 2 and 
ai t 2 control the mean variance, skewness (third central 
moment divided by the third power of the standard de- 
viation) and kurtosis of y s kow, which can be expressed 
in a closed form. Linearizing the dependence of these 
moments on the above parameters leads to a more com- 
plex system of equations. One is thus forced to obtain 
the parameters' values by multidimensional nonlinear fit- 
ting, losing the simplicity of the approximation (fTS|) . On 
the other hand, taking into account the skewness of the 
exact density of states allows to model better the exci- 
tation spectra with strongly non-homogeneous densities, 
such as E s — s 2 (square potential well) . This example is 
presented in Fig. [5] and compared with another skewness- 
sensitive approximation, the Edgeworth expansion (see 
e.g. [H|). The latter uses Chebyshev-Hermite polynomi- 
als to approximate a non-Gaussian probability distribu- 
tion, and thus does not guarantee the positivity of prob- 
ability density (Fig.O inset). In our comparison we have 
used the second-order Edgeworth expansion, which uses 
the same set of distribution moments as our fourth-order 
approximation. It is given by the formula 

_ e -v 2 n / #3(^3 H 4 (y) K4 H 6 (y)4 \ 
{V> V2^r \ 6a3 + 24a4 + 72^ J ' 

where H r is the r-th Chebyshev-Hermite polynomial, 
y = [M — Ki)/a, Kk is the fc-th cumulant of the ap- 
proximated distribution and a = ^/k2- The proposed 
skewness-sensitive approximation is as close, or closer 
(for low values of M) to the exact results as the Edge- 
worth expansion, and does not generate negative values. 
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Oi 




FIG. 5: Logarithms of exact and approximate values of 
u)(N, M, S) for N = 50 and 20 modes (classical particles, 
quadratic mode spacings). The inset shows the range of M 
values for which the Edgeworth expansion fails to preserve 
the positivity of the sum of states. 



IV. PHYSICAL APPLICATIONS 

In this section we discuss a number of physical prob- 
lems in which our exact iterative formulas and analytic 
approximations can be applied. 

The iterative expression @ can be used with arbitrar- 
ily discretized or real- valued mode excitations (measured 
or obtained numerically), such as isolated quantum dots 
with known number of electrons, often used to realize 
qubits (see e.g. [HI). The "modes" of a single dot are 
individual electronic configurations, and the mode ex- 
citations are the total energies of electrons inside the 
dot. A number of such dots can be thus treated as a 
system of classical independent particles, and its ther- 
modynamic properties described using the iterative for- 
mula © ■ Since it treats every system configuration sep- 
arately, this formula can also be used outside thermody- 
namics to calculate quantum properties which depend on 
the details of particular configurations, such as Hamilto- 
nian or thermal density matrix elements. 

The next problem to which our fourth-order approxi- 
mation can be applied is the entropy of a black hole and 
its relation to the area of the event horizon @, Q • It is 
well known that the appropriate ensemble in which to cal- 
culate the black hole entropy is the microcanonical one. 
This is true even for evaporating black holes, as they are 
not in the thermodynamic equilibrium [16j . Within the 
Loop Quantum Gravity theory, which quantizes geome- 
try, entropy as the function of the event horizon area will 
be proportional to the logarithm of the number of differ- 
ent microstates of quantized area observables which sum 
up to a given total area. As these observables are distin- 
guishable, we should use the classical particle statistics. 
The previous calculations of the microcanonical black 
hole entropy rely on the assumption that the excitation 
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FIG. 6: Comparison of the non-Gaussian approximations (our 
fourth-order and the second-order Edgeworth expansion) for 
the case of Fermi statistics and equal mode spacings, with 
N = 20 and 50 modes. 



spectrum of the area operator is evenly spaced 0, 
which is not true [l?} • Our results allow to extend the an- 
alytic approach to capture the full spectrum of the area 
operator excitations. Similarly, they can be applied to 
the statistics of photons in an isolated microwave cavity 
containing a number of distinguishable atoms [2j . 

Another problem for which our approximation is useful 
is the calculation of the nuclear state density i.e. the 
density p(M) of nuclear microstates which realize a par- 
ticular value M of the z component of the total angular 
momentum. Bethe approximated p(M) with a Gaussian 
distribution (using the Central Limit Theorem) and as- 
suming classical statistics for the particles constituting 
the nucleus This was a major drawback of his ap- 
proach, as the constituent particles of a nucleus are really 
fermions, and thus the Bethe approximation can be used 
only for small values of the total angular momentum of 
the nucleus [lsj . An improvement of this approximation 
uses the Edgeworth expansion (see e.g. 11]) matching 
the cumulants of p{M) calculated either in the canoni- 
cal or grand canonical ensemble @. As we have stated 
in Sec. MI CI the Edgeworth expansion does not guar- 
antee the positivity of probability density. Our method 
is free from this problem. Although both methods cap- 
ture higher moments of the density of states, the latter 
has the advantage that the resulting entropy fcs In Q is a 
polynomial function of the total excitation, and thus is 
easier to treat analytically. Additionally, the calculation 
of matched moments is not relegated to other statistical 
ensembles, but is done consistently in the microcanonical 
one — see (fT9|) and (j20|) . Consequently, as Fig. [6] shows, 
our approximation achieves a better fit to the exact cal- 
culations as compared with the second-order Edgeworth 
expansion (which utilizes the same set of moments as our 
approximation, but is nevertheless more complicated). 
The presented results are for Fermi statistics only, but 
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the comparison for the other two gives similar results. 

Finally, we note that the free energy F = E — TS = 
E — ksT In Q(M) becomes within our approximation a 
polynomial function of the total excitation M, which is 
particularly convenient when using the Ginzburg-Landau 
formalism (see e.g. |6|). Because it preserves the higher 
order M term in the exponent, our method can be par- 
ticularly useful when modeling systems containing small 
numbers of particles, in which the deviations from the 
Gaussian limit become stronger. It has this advantage 
also over the Edgeworth approximation, which (although 
taking into account higher moments) has only the M 2 
term in the exponent of the density of states, and there- 
fore does not lead to any free energy terms of higher than 
quadratic order. 

V. SUMMARY 

We have derived two exact expressions for the marginal 
or cumulative sum of states of classical, bosonic or 



fermionic particles, which are useful in practical calcu- 
lations. Exploiting a probabilistic analogy, we have then 
obtained an analytic fourth-order approximation to the 
normalized sum of states (density of states) u), which cap- 
tures its variance and kurtosis. The approximation can 
be extended to an analytic-numerical scheme which also 
matches the skewness of uj and is better suited to handle 
excitation spectra with strongly non-homogeneous den- 
sities. As shown numerically, for all three particle statis- 
tics our analytic formula is superior to the commonly 
used Gaussian one, which captures only the first and 
second moment of the density of states, and also to the 
Edgeworth expansion. Additionally, our approach em- 
ploys a simple exact method of calculating the moments 
of the microcanonical density of states, which requires 
less computational effort than the commonly used ap- 
proximations. To prove the above, we have tested our 
method numerically and discussed its applicability to var- 
ious physical systems. 
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